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Abstract 

Population  variability  and  uncertainty  are  important  features  of  biological  systems 
that  must  be  considered  when  developing  mathematical  models  for  these  systems.  In 
this  paper  we  present  probability-based  parameter  estimation  methods  that  account 
for  such  variability  and  uncertainty.  Theoretical  results  that  establish  well-posedness 
and  stability  for  these  methods  are  discussed.  A  probabilistic  parameter  estimation 
technique  is  then  applied  to  a  toxicokinetic  model  for  trichloroethylene  using  several 
types  of  simulated  data.  Comparison  with  results  obtained  using  a  standard,  deter¬ 
ministic  parameter  estimation  method  suggests  that  the  probabilistic  methods  are 
better  able  to  capture  population  variability  and  uncertainty  in  model  parameters. 
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1  Introduction 


Uncertainty  is  an  inherent  factor  in  mathematical  models  for  biological  sys¬ 
tems.  Model  equations  themselves  are  an  approximation  of  the  phenomena 
they  are  designed  to  model,  introducing  a  degree  of  uncertainty  that  is  dif¬ 
ficult  to  measure.  Further  simplifications  and  approximations  of  a  model  for 

*  Author  to  whom  correspondence  should  be  addressed. 

Email  addresses:  htbanks@eos.ncsu.edu  (H.T.  Banks), 
laura.k.potter@gsk.com  (Laura  K.  Potter). 

1  Current  address:  Scientific  Computing  and  Mathematical  Modeling,  Glaxo¬ 
SmithKline,  Research  Triangle  Park,  NC  27709  USA. 


Preprint  submitted  to  Mathematical  Biosciences 


13  September  2002 


Report  Documentation  Page 

Form  Approved 

OMB  No.  0704-0188 

Public  reporting  burden  for  the  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources,  gathering  and 
maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information, 
including  suggestions  for  reducing  this  burden,  to  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports,  1215  Jefferson  Davis  Highway,  Suite  1204,  Arlington 

VA  22202-4302.  Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person  shall  be  subject  to  a  penalty  for  failing  to  comply  with  a  collection  of  information  if  it 
does  not  display  a  currently  valid  OMB  control  number. 

1.  REPORT  DATE 

12  SEP  2002 

2.  REPORT  TYPE 

3.  DATES  COVERED 

00-00-2002  to  00-00-2002 

4.  TITLE  AND  SUBTITLE 

Probabilistic  methods  for  addressing  uncertainty  and  variablity  in 
biological  models:  Application  to  a  toxicokinetic  model 

5a.  CONTRACT  NUMBER 

5b.  GRANT  NUMBER 

5c.  PROGRAM  ELEMENT  NUMBER 

6.  AUTHOR(S) 

5d.  PROJECT  NUMBER 

5e.  TASK  NUMBER 

5f.  WORK  UNIT  NUMBER 

7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

North  Carolina  State  University, Center  for  Research  in  Scientific 
Computation, Raleigh, NC, 27695-8205 

8.  PERFORMING  ORGANIZATION 

REPORT  NUMBER 

9.  SPONSORING/MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

10.  SPONSOR/MONITOR'S  ACRONYM(S) 

11.  SPONSOR/MONITOR'S  REPORT 
NUMBER(S) 

12.  DISTRIBUTION/AVAILABILITY  STATEMENT 

Approved  for  public  release;  distribution  unlimited 

13.  SUPPLEMENTARY  NOTES 

The  original  document  contains  color  images. 

14.  ABSTRACT 

see  report 

15.  SUBJECT  TERMS 

16.  SECURITY  CLASSIFICATION  OF: 

17.  LIMITATION  OF 

ABSTRACT 

18.  NUMBER 

OF  PAGES 

41 

19a.  NAME  OF 

RESPONSIBLE  PERSON 

a.  REPORT 

unclassified 

b.  ABSTRACT 

unclassified 

c.  THIS  PAGE 

unclassified 

Standard  Form  298  (Rev.  8-98) 

Prescribed  by  ANSI  Std  Z39-18 


theoretical  and  computational  purposes  result  in  additional  layers  of  uncer¬ 
tainty.  Moreover,  many  biological  processes  are  subject  to  variability  that 
may  not  be  incorporated  into  a  mathematical  model.  Experimental  observa¬ 
tions  also  introduce  uncertainty  when  data  are  used  with  a  model  to  estimate 
parameters. 

Two  types  of  variability  that  are  common  in  biological  models  and  are  well- 
known  in  the  statistical  literature  are  intra-individual  and  inter-individual 
variability.  Intra-individual  variability  is  defined  as  variability  that  occurs 
within  a  given  individual  organism  or  biological  process.  This  type  of  variabil¬ 
ity  may  result  in  time-dependent  and/or  spatially-dependent  variation  within 
an  individual.  Biological  examples  of  such  variability  include  parameters  such 
as  body  weight,  blood  pressure,  fat  content  and  cell  membrane  permeabilities. 

A  second  type  of  variability  that  is  commonly  found  in  biological  modeling  is 
inter-individual  variability.  This  type  of  variability  results  from  variations  in 
individuals  across  a  population.  Biological  models  that  are  based  on  behavior 
or  phenomena  over  a  population  are  almost  always  subject  to  inter-individual 
variability.  This  is  especially  the  case  when  a  model  is  designed  to  predict  or 
explain  experimental  observations  that  are  collected  from  multiple  individuals. 

It  is  reasonable  to  expect  that  different  individuals  of  a  population  would 
possess  different  values  for  biological,  physical  and  chemical  parameters.  These 
parameters  would  then  take  on  a  range  of  values  across  the  population,  so  that 
each  parameter  would  be  associated  with  a  probability  distribution  that  would 
mathematically  describe  this  variation.  Using  data  from  multiple  individuals, 
the  resulting  probability  distributions  can  be  estimated  with  inverse  problem 
techniques. 

Examples  of  biological  parameters  that  are  often  subject  to  inter-individual 
variability  include  growth  and  death  rates,  susceptibility  to  infection,  efficacy 
of  vaccines  and  other  prophylactics,  and  age.  Note  that  each  of  the  examples 
given  above  for  intra-individual  variability  may  also  involve  inter-individual 
variability  depending  on  the  type  of  model  and  experimental  observations. 
Similarly,  each  of  the  examples  for  inter-individual  variability  also  may  be 
subject  to  intra- individual  variability. 

The  motivating  example  we  consider  here  is  a  toxicokinetic  model  for  the 
systemic  transport  of  the  environmental  contaminant  trichloroethylene  (TCE). 
TCE  is  a  solvent  that  has  been  used  widely  in  industry  as  a  metal  degreasing 
agent,  and  is  now  a  common  soil  and  groundwater  contaminant.  This  highly 
fat-soluble  compound  is  rapidly  absorbed  into  the  bloodstream,  and  has  been 
shown  to  accumulate  in  the  adipose  (fat)  tissue  of  humans  and  animals  [1,2]. 
Known  and  suspected  toxic  effects  of  TCE  and  its  metabolites  in  laboratory 
animals  and/or  humans  include  acute  effects  such  as  dizziness,  drowsiness, 
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headaches  and  fatigue,  as  well  as  chronic  effects  such  as  developmental  defects 
and  lung,  kidney  and  liver  tumors  [3-10]. 

Toxicokinetic  models  are  used  in  the  overall  risk  assessment  process  for  toxic 
compounds  to  help  quantify  the  expected  risk  of  toxicity  to  humans  as  a 
function  of  the  level  of  exposure  to  the  given  chemical.  In  particular,  phys¬ 
iologically  based  pharmacokinetic  (PBPK)  models  predict  the  effective  dose 
level  of  a  toxic  compound  that  is  delivered  to  the  “target”  tissues  (i.e.,  tis¬ 
sues  that  experience  toxic  effects)  for  a  given  external  exposure  level.  PBPK 
models  are  compartmental  models  that  describe  the  systemic  transport  of 
a  compound  through  the  tissues  and  organs,  including  the  dynamics  of  up¬ 
take,  tissue  distribution,  metabolism  and  elimination.  The  resulting  model  is 
a  system  of  ordinary,  partial  and/or  delay  differential  equations,  with  each 
equation  representing  the  dynamics  of  tissue  concentrations  in  a  particular 
tissue  or  organ. 

Several  PBPK  models  have  been  developed  for  TCE  (e.g.,  see  [11-15]).  A 
majority  of  these  models  make  use  of  the  standard  “perfusion-limited”  com¬ 
partmental  model  for  each  of  the  modeled  tissues  and  organs  (see  Section  4 
for  a  description  of  the  perfusion-limited  model).  In  [16]  and  [17],  three  PBPK 
models  for  TCE  are  developed  and  compared,  each  with  a  different  submodel 
for  the  adipose  tissue  compartment. 

As  discussed  in  [16],  preliminary  simulations  indicated  that  a  perfusion- limited 
adipose  tissue  compartment  does  not  appear  to  sufficiently  capture  the  dynam¬ 
ics  of  TCE  accumulation  in  fat  as  seen  in  experimental  data.  Moreover,  adipose 
tissue  is  known  to  have  highly  heterogeneous  physiological  properties,  includ¬ 
ing  significant  variations  in  fat  cell  size,  lipid  distribution,  blood  flow  rates  and 
cell  membrane  permeabilities  [18-21],  These  characteristics  further  suggest 
that  the  “well-mixed,”  rapid  equilibrium  assumptions  of  the  perfusion-limited 
model  may  be  inappropriate  for  describing  the  disposition  of  fat-accumulating 
compounds  such  as  TCE  in  adipose  tissue. 

To  better  capture  the  dynamics  of  TCE  in  fat  tissue,  a  spatially  varying  axial 
dispersion  model  was  developed  [16]  to  address  the  intra-individual  variabil¬ 
ity  that  results  from  the  heterogeneous  lipid  distribution  and  physiological 
characteristics  of  adipose  tissue.  This  variability  is  built  into  the  adipose  com¬ 
partmental  model  with  a  special  axial  dispersion  term,  where  the  “dispersion” 
coefficient  is  a  measure  of  the  degree  of  intra- individual  variability  that  occurs 
in  the  fat. 

In  addition  to  the  intra-individual  variability  that  appears  to  affect  TCE  con¬ 
centrations  in  fat  tissue,  inter-individual  variability  also  plays  a  major  role 
in  toxicokinetic  models  in  general.  Current  technology  almost  always  requires 
that  measurements  of  chemical  concentrations  in  tissues  over  time  must  be 
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taken  from  multiple  individuals.  This  experimental  necessity  immediately  in¬ 
troduces  inter-individual  variability  into  the  measured  observations,  and  must 
be  considered  in  the  development  of  mathematical  models. 


As  biological  models  have  become  more  widely  utilized  and  influential  in  a 
variety  of  fields,  the  need  to  account  for  variability  and  uncertainty  in  mod¬ 
eling  has  been  recognized.  Markov  Chain-Monte  Carlo  methods  have  been 
developed  to  address  issues  of  variability  and  uncertainty,  and  these  methods 
have  been  applied  to  PBPK  models  as  a  part  of  the  parameter  estimation  pro¬ 
cess.  Monte  Carlo  methods  are  based  on  a  Bayesian  statistical  approach  that 
involves  the  use  of  experimental  data  to  update  estimates  of  a  hypothesized 
“prior”  probability  distribution  for  one  or  more  model  parameters.  Examples 
of  Monte  Carlo  methods  applied  to  PBPK  models  can  be  found  in  [22-31]. 


An  alternative,  probability-based  method  has  been  developed  to  incorporate 
uncertainty  and  variability  in  mathematical  models.  This  method,  which  is 
discussed  in  [32-34]  and  is  detailed  in  Section  2,  is  centered  around  a  proba¬ 
bilistic  parameter  estimation  approach  that  involves  the  estimation  of  proba¬ 
bility  distributions  for  model  parameters.  Well-known  theoretical  results  from 
probability  theory  establish  the  theoretical  soundness  of  this  technique,  which 
can  be  implemented  computationally  in  a  straightforward  manner. 


A  distinct  advantage  of  this  approach  over  the  Monte  Carlo-based  methods  is 
an  added  level  of  flexibility  in  choosing  the  prior  probability  distributions.  As 
we  discuss  in  Section  2,  these  probability-based  methods  can  be  used  with  pre¬ 
selected  prior  distributions  as  with  Monte  Carlo  methods,  or  they  may  be  used 
with  weighted  sums  of  Dirac  delta  measures  that  do  not  assume  a  fixed  form 
for  the  probability  distribution  functions.  A  version  of  this  method  has  been 
applied  to  a  population  model  for  mosquitofish  in  rice  paddies,  and  was  used 
to  successfully  describe  fish  population  dynamics  by  estimating  distributed 
growth  rate  functions  using  aggregate  experimental  data  [35]. 


In  this  paper  we  present  probability-based  parameter  estimation  methods  for 
incorporating  uncertainty  and  variability  into  biological  models.  These  meth¬ 
ods  are  general  and  may  be  applied  to  a  wide  variety  of  models  to  account 
for  various  types  of  model  uncertainty  as  we  have  outlined  here.  In  Section  2 
we  formulate  these  probabilistic  parameter  estimation  methods  in  a  general 
setting.  We  address  related  theoretical  issues  in  Section  3,  establishing  the 
wcll-posedness  of  the  resulting  parameter  estimation  process.  Finally,  imple¬ 
mentation  of  the  methods  in  the  context  of  a  toxicokinetic  model  for  TCE  is 
discussed  in  Section  4. 
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2  General  problem  formulation 


Suppose  that  a  biological  process  is  described  by  the  parameter-dependent 
system 


(1) 


where  y  represents  the  state  of  system,  q  is  the  vector  of  parameters  on  which 
the  state  depends,  and  /  represents  the  dynamics  of  the  system  in  the  form 
of  ordinary,  partial  or  functional  differential  equations.  Experimental  data 
z  =  {zi},  i  =  1  are  given  that  correspond  to  complete  or  partial 

observations  Oy(U ;  q)  of  the  state. 

The  model  parameters  q  are  estimated  in  a  deterministic,  least  squares  setting 
by  minimizing  the  objective  function 

Nt 

J{q,z)  =  '52\Oy(ti]q)-  Zil2  (2) 

i= 1 


over  q  G  Q  subject  to  (1),  where  Q  is  the  space  of  admissible  parameters  and 
O  is  the  observation  operator. 

In  standard  least  squares  estimation  problems,  the  space  Q  is  usually  de¬ 
fined  as  a  compact  subset  of  Mn  for  some  positive  integer  n,  so  that  q  = 
(<7i,  q-2,  ■  ■  • ,  qn)  is  a  vector  with  real  components.  This  assumption  requires 
each  parameter  to  be  a  constant,  which  may  not  be  reasonable  for  parameters 
and  experimental  data  that  are  subject  to  inter-individual  variability.  Indeed, 
if  experimental  observations  are  collected  from  multiple  individuals  in  a  pop¬ 
ulation,  then  one  must  think  of  the  constant  parameters  q  G  Mn  as  average 
values  over  the  sampled  population.  This  approximation  may  be  appropriate 
for  some  parameters  that  do  not  vary  to  a  large  extent  across  individuals, 
but  in  many  cases  these  “mean”  value  approximations  may  lead  to  inaccurate 
parameter  estimates  and  an  inaccurate  description  of  the  population.  This  is 
especially  true  in  situations  when  subpopulations  are  described  by  different 
parameter  values,  or  means,  variances,  etc. 

Such  population-dependent  variability  in  model  parameters  can  be  incorpo¬ 
rated  into  least  squares  estimation  problems  using  a  probability-based  formu¬ 
lation.  We  assume  that  the  model  parameters  q  are  realizations  of  random 
variables  with  probability  distributions  P  that  vary  over  the  population,  so 
that  P  belongs  to  a  probability  space  Q  that  may  be  infinite  dimensional.  As 
in  [33],  we  define  the  set  V(Q)  of  all  probability  distributions  on  the  admissible 


5 


parameter  space  Q  and  seek  a  probability  distribution  function  PeQC  P(Q) 
that  minimizes  the  objective  function 


Nt 


J(P,z)  =  'Z\£[Oy(ti;q)\P]-zi\ 


i=  1 


(3) 


over  Q  C  V(Q),  where  £j,  i  —  1, . . . ,  Nt  are  observations  corresponding  to  the 
expected  value 

£[Oy(t,;q)\P-}  =  J  Oy(U,q)dP-(q)  (4) 

Q 


for  some  P*  G  Q  C  V(Q).  For  simplicity  we  often  choose  Q  =  V(Q),  but 
this  is  not  essential  and  one  may  readily  restrict  the  family  of  admissible 
distributions  in  certain  formulations. 

Depending  on  the  choice  of  the  set  Q  C  V(Q)  of  probability  distributions, 
this  method  may  be  implemented  with  pre-determined  “prior”  probability 
distributions  (as  with  the  Monte  Carlo  method),  or  it  may  be  used  without 
the  pre-specification  of  a  particular  probability  distribution.  For  the  case  when 
there  is  a  reasonable  expectation  that  a  parameter  varies  across  the  population 
in  a  manner  similar  to  a  given  probability  distribution,  the  set  Q  can  be  chosen 
as  the  space  of  those  distribution  functions  (e.g.,  log  normal  distributions) 
defined  over  the  admissible  parameter  space  Q.  For  this  type  of  formulation, 
the  distribution  functions  are  uniquely  determined  by  their  parameterizations 
q  (e.g.,  q  =  (/x,  cr) ,  the  mean  and  standard  deviation),  and  hence  may  be 
estimated  by  minimizing 

Nt 

J(q,z)  =  \£[°y(<ti-,q)\P(q)]  -  Zi\2  (5) 

i=  1 


over  the  space  Q  of  admissible  parameterizations  q,  where  Q  is  given  as  the 
set  of  probability  distributions  of  the  pre-specihed  form  with  the  parameteri¬ 
zations  q. 

For  example,  if  it  is  believed  that  a  parameter  can  be  approximated  by  a 
log  normal  distribution,  then  Q  is  given  as  the  set  of  all  log  normal  distribu¬ 
tions  defined  over  Q.  For  this  particular  formulation,  the  estimation  problem 
is  solved  by  minimizing  (5)  over  the  space  Q  of  admissible  mean  and  stan¬ 
dard  deviation  parameters  q  =  (n,  a).  This  approach  has  been  implemented 
in  [36]  where  Q  is  defined  as  a  set  of  bitruncated  normal  distributions  with 
certain  specified  properties,  and  the  estimated  parameters  q  are  the  mean  and 
a  modified  standard  deviation  (see  [36]  for  details). 
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If  it  is  not  possible  to  predict  the  expected  form  of  the  probability  distributions 
a  priori ,  this  method  also  may  be  used  without  the  specification  of  prior  distri¬ 
butions.  In  this  case,  Q  =  V(Q)  may  be  chosen  as  the  space  of  all  probability 
distributions  defined  on  Q.  For  computational  purposes,  the  estimation  prob¬ 
lem  may  then  be  implemented  using  finite  dimensional  approximations  to  the 
original  infinite  dimensional  problem.  First  we  define  the  infinite  dimensional 
set 


Vo (Q)  =  {Pe  V(Q)  :  P  =  'ZPi^ke  N+,qj  e  Q0,Pj  >  0 =  l},(6) 


3= 1 


3  =  1 


where  Q0  =  {qj}'^=1  is  a  given  countable,  dense  subset  of  the  parameter  space 
Q  and  5qj  is  the  Dirac  delta  distribution  with  mass  at  qj  G  Q.  In  other  words, 
Vo  (Q)  is  the  set  of  probability  distributions  on  Q  that  have  finite  support  in 
Q0.  We  then  define  the  finite  dimensional  set 


M 

VM  =  {PeV0(Q):P  =  Y.ViK) 

3= 0 

which  we  use  to  define  a  family  of  finite  dimensional  approximation  problems. 
That  is,  for  fixed  {q0,  q2, . . . ,  ?m}  in  Q0  with  PM  =  T,jLoPj8qj  e  VM ,  we 
minimize  the  objective  function 


Nt 


J(PM,z)  =  \S[Oy(ti]  q)\PM]  -  Zi 


i— 1 


Nt 

=E 


i= 1 


E  °y(tt,  Qj  )pj  -  Zi 

3=0 


(7) 

(8) 


over  the  set  VAI .  These  precise  definitions  are  necessary  to  obtain  a  well- 
posed  estimation  problem,  as  we  discuss  in  Section  3.  Note  that  the  problem 
of  minimizing  the  objective  function  (8)  corresponds  to  solving  a  constrained 
quadratic  programming  problem  for  (p0, . . .  ,pM)  with  the  constraints  Pj  >  0, 
Y^jLoPj  =  1  (see  Section  4.3).  There  currently  exist  a  number  of  acceptable 
computational  methods  to  solve  such  problems  which  are  again  special  cases 
of  choosing  an  a  priori  parameterization  set  (Q  =  VAI  in  this  case)  and 
optimizing  over  admissible  parameterizations 


M 

Q  =  {q=  (P(h  . . .  ,pM)  :  Pj  >  0 ,Y,Pj  =  !}■ 

3=0 
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3  Theoretical  results 


In  order  to  address  theoretical  issues  related  to  the  inverse  problems  discussed 
in  Section  2,  we  need  to  define  a  suitable  metric  for  the  probability  spaces 
V{Q).  Using  the  Prohorov  metric  and  well-known  results  from  probability 
theory,  we  establish  a  theoretical  framework  that  allows  us  to  prove  method 
stability  for  our  probability-based  parameter  estimation  problems. 


3.1  The  Prohorov  metric 


As  in  [37],  the  Prohorov  metric  is  defined  on  the  space  of  probability  measures 
V{Q)  on  the  Borel  subsets  of  Q,  where  Q  is  a  complete  metric  space  with 
metric  d.  The  Prohorov  metric  p  :  V{Q)  x  V(Q)  — >  M+  is  defined  by 

p(Pi,  P2)  =  inf{e  >  0  :  P1[F]  <  P2[Fe }  +  e;  F  closed;  F  C  Q}, 


where 

F£  =  {q  E  Q  :  d(q ,  q)  <  e;  q  G  F}. 

It  is  well  known  that  p  is  a  metric  on  the  space  V(Q),  and  that  this  metric 
space  ( V(Q),p )  is  complete.  Moreover,  (' P(Q),p )  is  compact  for  all  compact 
sets  Q. 

Another  well-known  result  [37]  establishes  equivalent  conditions  for  conver¬ 
gence  of  probability  distributions  in  the  Prohorov  metric.  Assuming  that  (Q,  d) 
is  complete,  then  the  following  convergence  statements  for  Pk,P  €  V(Q)  are 
equivalent: 

(i)  p(Pk,  P )  -*•  0. 

(ii)  JQf(q)dPk{q)  ->  SQ  f(q)dP(q)  for  all  bounded  and  uniformly  continuous 
functions  /  :  Q  >  M. 

iii)  Pk[A)  — >  P[A]  for  all  Borel  sets  A  C  Q  with  P[dA]  =  0,  where  d A  is  the 
boundary  of  A. 


3.2  Stability  of  the  general  parameter  estimation  problem 


Banks  and  Bihari  [33]  have  addressed  theoretical  issues  related  to  probability- 
based  estimation  problems.  Using  the  Prohorov  metric,  they  studied  the  con¬ 
vergence  properties  of  sequences  of  probability  distributions  in  V(Q).  These 
results  were  then  applied  to  a  sequence  of  minimizers  for  finite  dimensional 
approximations  to  the  estimation  problem  for  (3).  Here  we  summarize  their 
findings  as  they  relate  to  the  inverse  problems  described  in  Section  2. 


As  discussed  in  [33],  it  follows  that  if  the  mapping  q  — >  Oy(U ;  q)  is  continuous, 
then  the  convergence  p(Pk,  P)  — >  0  in  the  Prohorov  metric  is  equivalent  to 

S[Oy(ti]  q)\Pk\  ->  £[eh/(U;g)|-P], 


and  hence  the  map 

Nt 

P^J(P)  =  Y,\£[Oy(ti,q)\P]-zi\2 

i= 1 

is  continuous  in  the  p  topology.  Moreover,  if  the  space  Q  is  compact  we  have 
that  (V(Q),p)  is  a  compact  metric  space,  which  along  with  the  continuity  of 
the  map  P  J(P)  guarantees  the  existence  of  a  minimizer  over  V{Q)  for  the 
estimation  problem  associated  with 

Nt 

mmJlT,  i)  =  ]T  \£[Oy(ti ;  q)\P]  -  z^2  .  (9) 

r(E.  r{Q)  i=\ 

In  addition  to  establishing  the  existence  of  a  solution  for  the  inverse  prob¬ 
lem  (9),  Banks  and  Bihari  [33]  developed  results  related  to  method  stability  for 
this  problem.  Using  finite  dimensional  approximation  techniques,  they  show 
in  Theorem  4.1  that  the  solutions  for  (9)  depend  continuously  on  the  data 
(see  [33]  for  a  complete  discussion).  Moreover,  any  sequence  of  minimizers  of 
the  finite  dimensional  problems  for  (7)  converge  in  the  Prohorov  metric  to 
a  minimizer  for  the  original  infinite  dimensional  problem  (9).  This  theorem 
makes  use  of  the  result  they  prove  in  Theorem  3.1  [33]  that  the  set  Vq{Q)  as 
in  (6)  is  dense  the  space  V(Q)  with  respect  to  the  Prohorov  metric  p. 

In  demonstrating  the  convergence  of  solutions  for  the  family  of  finite  dimen¬ 
sional  problems  (7),  the  result  established  in  Theorem  4.1  of  [33]  also  provides 
a  computational  framework  for  solving  the  general  parameter  estimation  prob¬ 
lem  (9)  without  specifying  prior  probability  distributions.  Using  discrete  Dirac 
delta  measures,  for  sufficiently  large  M  we  may  approximate 

f  r  M  M 

/  Oy{t ;  q)dP(q)  ~  (i)P3dqM  (q)dq  =  ^  Oy(t ;  qf)pj: 

Q  Q  J=°  '  3=° 

which  then  allows  us  to  approximate  the  infinite  dimensional  inverse  prob¬ 
lem  (9)  by  the  finite  dimensional  approximation  (7). 


4  Application  to  a  toxicokinetic  model  for  TCE 


In  this  section  we  apply  the  methods  developed  in  Section  2  to  a  toxicoki¬ 
netic  model  for  trichloroethylene.  We  utilize  the  probability-based  estimation 
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technique  involving  (3)  and  we  compare  the  results  to  those  obtained  with  a 
traditional  deterministic  method  for  (2). 

Here  we  develop  and  test  several  estimation  problems  for  the  TCE  model  with 
simulated  data  that  qualitatively  and  quantitatively  match  the  experimental 
data  in  [38],  and  we  demonstrate  which  strategies  and  observations  are  most 
successful  at  capturing  and  predicting  the  dynamics  of  TCE  in  adipose  tis¬ 
sue.  This  in-depth  study  of  the  parameter  estimation  process  with  simulated 
data  is  important  for  testing  and  understanding  the  capability  of  these  esti¬ 
mation  techniques ,  and  is  a  necessary  first  step  before  use  of  the  methods 
with  experimental  data  containing  additional  and  generally  unknown  sources 
of  variability. 

The  results  that  we  present  here  demonstrate  clear  advantages  for  the  probability- 
based  method  when  the  experimental  observations  are  subject  to  inter-individual 
variability.  In  addition,  our  studies  of  the  parameter  estimation  problem  for 
the  TCE  model  illustrate  ways  in  which  both  the  quantity  and  the  quality 
of  experimental  data  have  a  major  impact  on  the  effectiveness  of  parameter 
estimation  techniques. 


4-1  Overview  of  the  TCE  model 


Here  we  provide  an  overview  of  the  PBPK-hybrid  model  for  TCE  as  developed 
in  [16,36].  This  model  utilizes  standard  physiologically  based  pharmacokinetic 
compartmental  equations  for  various  non-fat  tissues.  The  fat  tissue  compart¬ 
ment  is  described  with  a  spatially  varying  dispersion  model,  and  is  designed 
specifically  to  capture  the  intra-individual  variability  that  results  from  the 
heterogeneous  physiology  of  fat. 

The  most  commonly  used  compartmental  model  in  PBPK  modeling  is  the 
perfusion-limited,  or  flow-limited  compartment.  This  model  is  based  on  sim¬ 
ple  mass  balance  principles  and  assumptions  of  rapid  equilibrium  and  spatial 
uniformity.  Moreover,  it  is  assumed  that  the  blood  flow  rate  to  the  tissue 
is  much  slower  than  the  rate  of  transport  of  the  compound  across  cell  mem¬ 
branes.  The  resulting  equation  for  the  tissue  concentration  C  of  the  compound 
is  given  by 

yd^  =  QbKCin(t)  ~  Cout(t)), 

where  V  is  the  volume  of  the  tissue,  Qu  is  the  volumetric  blood  flow  rate 
to  the  tissue,  and  Cin  and  Cout  are  the  concentrations  of  compound  entering 
and  leaving  the  tissue  respectively  (see  [39]).  Under  standard  assumptions,  the 
concentration  Cout  is  equal  to  the  concentration  C  of  compound  in  the  tissue 
divided  by  the  blood:tissue  partition  coefficient  [39]. 
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PBPK  Model  for  Inhaled  TCE 


Fig.  1.  Schematic  of  PBPK  model  for  inhaled  TCE  in  Long-Evans  rats. 

For  many  tissues  and  compounds  of  interest  the  perfusion-limited  compart- 
mental  model  is  adequate  to  describe  the  dynamics  of  such  compounds  inside 
the  tissues.  In  the  case  of  highly  lipophilic  substances  such  as  TCE,  however, 
the  standard  models  may  not  accurately  capture  the  transport  of  these  chem¬ 
icals  in  the  adipose  tissue.  As  discussed  in  Section  1,  the  highly  heterogeneous 
physiology  of  fat  tissue  appears  to  have  a  major  influence  on  the  behavior 
of  TCE  in  fat.  Using  a  PBPK  model  for  TCE  in  Long-Evans  rats  with  a 
perfusion-limited  fat  compartment  [38]  (see  Figure  1  for  a  model  schematic), 
model  simulations  suggested  that  the  standard  model  indeed  does  not  capture 
the  concentration  profile  of  TCE  in  adipose  tissue  as  seen  in  experimental 
data  [16]. 

To  account  for  the  spatial  variation  in  TCE  fat  concentrations  as  suggested 
by  the  physiology  of  adipose  tissue,  an  axial  dispersion  model  was  developed 
to  replace  the  perfusion-limited  fat  tissue  compartment.  This  model  is  based 
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Fig.  2.  Geometric  representation  of  an  adipocyte-capillary  unit  in  adipose  tissue. 
The  adipocyte  region  (A)  is  represented  by  a  sphere  and  is  immersed  in  the  inter¬ 
stitial  fluid  (I).  The  capillary  or  blood  region  (B)  is  a  cylindrical  tube  that  wraps 
around  the  adipocyte.  Coordinates  are  in  spherical  coordinates  (r,  0,  </>). 


directly  on  the  structure  of  fat  tissue,  which  consists  primarily  of  spherical, 
lipid-containing  cells  called  adipocytes.  Each  adipocyte  is  in  contact  with  one 
or  more  capillaries  [40]  and  is  immersed  in  interstitial  fluid.  Figure  2  depicts 
the  geometric  representation  of  an  adipocyte-capillary  unit  as  used  in  the  dis¬ 
persion  model.  There  are  three  subcompartments  in  the  model  that  represent 
the  adipocyte  region  (A),  the  capillary  or  blood  region  (B)  and  the  interstitial 
fluid  (I). 

The  model  equations  are  based  on  an  axial  dispersion  model  developed  by 
Roberts  and  Rowland  [41]  for  the  liver.  A  key  feature  of  their  model  is  its 
aggregate  structure ,  using  a  single  cellular  unit  with  the  dispersion  term  to 
represent  the  intra-individual  variability  that  occurs  across  the  millions  of 
cells  in  the  tissue.  As  detailed  in  [16],  we  have  adapted  their  model  to  describe 
the  geometry  of  adipose  tissue  and  the  transport  of  TCE  within  the  fat.  The 
resulting  system  of  partial  differential  equations  is  given  by 


y  dC'B  _  Vb  d 
dt  r  2  sin  0  dcj) 

+  \iHBi{fiCi(90)  —  / bCb ) 

+  ^a^baU aCa{9  o)  —  JbCb) 


sm< 


~Db  dCB 

r  2  dcj> 


-vC 


B 


(10) 
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+ vCs{t’  v  = m (11) 

dCj_  VjDj  r  i  d2Cr  id  f .  dCjY 

1  dt  r'i  sin2  0  d92  sin  <p  d<f>  \  d<p  ) 

+  ^o($)Xb(0M/Rb/(/bC,b  —  fiCi ) 


+  Via^/aCa  —  fiCi)  (13) 

CI(t,9,(j>)  =  CI(t,9  +  2n,(j))  (14) 

BC1  BC 

-Bl(t,9,<l>)  =  -^(t,9  +  2tt,0)  (15) 

Cj(t,  9, 0)  <  oo  (16) 

Ci(t,9, 7r)  <  oo  (17) 


dCA  VaDa  \  1  &CA  Id  (  ,  dCA\ 
dt  Tq  sin2  <f>  d 92  sin  <f>  d(p  y  d<p  ) 
+  ^e0(9)xB(4>)^A^BA(fBCB  ~  f aCa ) 


+  ViaUiCi  —  fACA)  (18) 

CA{t,9,<j>)  =  CA{t,9  +  27t»  (19) 

BC1  BC 

-^(t,9,<j))  =  ^(t19  +  2n,(j))  (20) 

CU(M,0)<oo  (21) 

CA(t,9, 7r)  <  oo.  (22) 


The  capillary  equation  (10)  describes  the  transport  of  TCE  in  the  capillary 
region  of  the  adipose  tissue  and  utilizes  the  dispersion  term 

Vr  d  .  T>b  dC'B 

- — 7  7T7  W - W~ 

r2sm<pd(p  [  r'2  dcp 

with  dispersion  coefficient  Tffi.  This  term  accounts  for  the  variability  in  phys¬ 
iological  properties  that  occurs  across  the  population  of  fat  cells,  with  a  large 
dispersion  coefficient  indicating  a  high  degree  of  variability.  Mathematically, 
the  dispersion  term  is  equivalent  to  a  standard  diffusion  term,  although  the 
dispersion  term  is  used  specifically  to  approximate  the  observed  physiologi¬ 
cal  phenomena  of  varying  path  lengths,  flow  velocities  and  compound  transit 
times  that  occur  within  a  tissue. 

The  boundary  conditions  (11)  and  (12)  connect  the  adipose  capillary  region 
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to  the  systemic  arterial  and  venous  blood  compartments  using  flux  balance. 
Transport  of  TCE  between  the  capillary  region  and  the  other  two  adipose 
subcompartments  (interstitial  and  adipocyte)  is  modeled  in  the  PDE  (10). 
The  variables  Cb (t),  Cj{t )  and  CaX)  denote  concentrations  of  TCE  in  the 
capillary,  interstitial  and  adipocyte  regions  respectively,  while  Ca(t )  and  Cv(t) 
represent  the  systemic  arterial  and  venous  blood  concentrations  of  TCE. 

The  interstitial  region  is  modeled  with  the  two-dimensional  PDE  (13)  and 
boundary  conditions  (14)  -  (17).  The  adipocyte  region  equations  (18)  -  (22) 
are  similar  in  structure  to  the  interstitial  equations,  and  describe  the  diffusion 
of  TCE  around  the  surface  of  the  adipocyte  as  well  as  the  transport  of  TCE 
between  the  three  adipose  subcompartments.  The  boundary  conditions  are 
standard  periodic  and  finiteness  boundary  conditions  that  are  commonly  used 
for  the  diffusion  equation  on  a  spherical  domain.  A  detailed  derivation  and 
description  of  the  dispersion  model  is  given  in  [16,36]. 

Adipose  model  parameters  include  the  dispersion  coefficient  T>b  (m2/hour); 
diffusion  coefficients  Dj  and  Da  (m2/hour);  the  fractions  fs,  fi ,  f  a  of  unbound 
TCE  in  each  adipose  region;  cell  membrane  permeability  coefficients  /pba,  l1 ia , 
Ubi  (liters /hour);  blood  flow  parameters  v  (m/hour)  and  T\  and  inter-region 
transport  parameters  A/  and  A  a- 

The  adipose  model  equations  (10)  -  (22)  are  coupled  with  standard  com- 
partmental  equations  for  the  lung,  arterial  blood,  venous  blood,  liver,  brain, 
kidney,  muscle  and  remaining  non-fat  tissue  to  obtain  a  whole-body  PBPK- 
hybrid  model.  Uptake  of  TCE  is  via  inhalation  into  the  lungs,  and  metabolism 
is  modeled  with  a  Michaelis-Menten  term  in  the  liver.  The  resulting  equations 
are  given  by 


—  QmPm/ Pm  +  QiPt/ Pt  +  Q fC b(',  ^  ~  £2)  +  QbrCbr/ Pbr 


+  QiQ/ Pi  +  QkCk/ Pk  —  QcCv 

QcCv  +  QpCc 

a~  ^  +  t 

(23) 

(24) 

Vmd^  =  Qm{Ca-Cm/P,n) 

(25) 

(1C 

v,dt‘=Qt(c„~c,iPt) 

(26) 

Vbr  ^  —  Qbripa  Cbr/Pbr ) 

(27) 

,  VmaxCi/ Pi 

V'dt=Ql(C‘  c,/ P,)  kM+C,/P, 

(28) 

VkdCJ=Qt(Ca  Ck/P„), 

at 

(29) 
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where  Cv(t),  Cbr(t),  Ck(t),  Cm(t),  Ci(t )  and  Ct(t)  denote  TCE  concentrations 
in  the  venous  blood,  brain,  kidney,  muscle,  liver  and  remaining  tissue  com¬ 
partments,  respectively.  The  chamber  air  concentration  Cc(t)  is  specified  as 
part  of  the  experiment  and  is  used  as  a  forcing  function  in  the  arterial  blood 
equation  (24).  For  the  results  we  present  in  this  paper,  we  set  the  chamber  air 
concentration  to  2000  parts  per  million  TCE  for  one  hour,  followed  by  zero 
ppm  TCE  until  the  final  time  tf  (in  hours). 

Model  parameters  include  tissue  volumes  V  (in  liters),  volumetric  blood  flow 
rates  to  the  tissues  Q  (liters/hour),  and  blood:tissue  partition  coefficients  P, 
each  labeled  with  a  subscript  corresponding  to  the  appropriate  tissue.  The  car¬ 
diac  output  and  ventilation  rates  (in  liters/hour)  are  denoted  by  Qc  and  Qb 
respectively,  and  the  blood:air  partition  coefficient  is  denoted  as  Pb .  The  stan¬ 
dard  Michaelis-Menten  metabolic  parameters  are  denoted  by  vmax  (mg/hour) 
and  h'M  (mg/liter).  See  [16,36]  for  complete  discussion  of  the  model  equations 
and  parameters. 

Theoretical  results  relating  to  well-posedness  of  the  whole-body  PBPK-hybrid 
model  are  presented  in  [42],  In  particular,  we  have  shown  the  existence  of  a 
unique  weak  solution  for  a  general  class  of  nonlinear  parabolic  equations  that 
includes  the  TCE  model  as  a  special  case.  Moreover,  we  established  the  well- 
posedness  of  the  deterministic  estimation  problem  for  the  TCE  model,  and 
in  [36]  we  have  addressed  the  well-posedness  of  probability-based  parameter 
estimation  methods  applied  to  the  TCE  model.  Numerical  methods  and  sim¬ 
ulations  for  this  model  with  deterministic  parameters  are  given  in  [17],  and 
results  for  the  standard  PBPK  models  are  compared  to  those  for  the  PBPK- 
hybrid  model. 


4-2  Deterministic  parameter  estimation  methods  for  the  TCE  model 


In  this  section  we  present  results  for  the  standard  deterministic  parameter 
estimation  problem 


Nt 

min  J(q,  z)  =  Q)  ~  zi 

qeQ  i= i 


(30) 


applied  to  the  TCE  model,  where 

y(t)  =  [< CB(t ),  C7(t),  CA{t)i  Cv(t),  Cbr{t ),  Ck(t),  Cm(t),  Ct(t),  G(t)]T, 

q  denotes  the  vector  of  unknown  parameters  in  the  admissible  parameter  space 
Q,  the  observations  are  denoted  by  z^,  i  =  1, . . . ,  Nt,  and  O  is  the  observation 
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operator.  Here  the  variable  y(t)  is  subject  to 

=  /(*>J/(V);?), 

which  we  use  as  abbreviated  notation  for  the  TCE  PBPK-hybrid  model  (10)  - 
(29). 

4-2.1  Simulated  observations 

The  observations  Zi,  i  —  1, . . . ,  Nt  and  the  observation  operator  O  are  defined 
so  that  they  correspond  to  the  types  of  experimental  observations  used  in 
the  experiments  conducted  by  Evans  et  ah  [38].  The  data  that  they  collected 
include  measurements  of  TCE  concentrations  in  homogenized  samples  of  fat 
tissue.  Therefore  we  define  the  observation  operator 


Oy(ti\  q)  =  CA(ti',  q)  (31) 

where  CUfe?)  is  the  mean  concentration  of  TCE  over  the  adipocyte  region 
(see  [17]  for  a  precise  definition  and  finite-dimensional  approximation).  Simi¬ 
larly,  we  define  the  simulated  observations 

Zi  =  CA{ti,  q*)  (32) 

for  i  —  1, . . . ,  Nt  and  for  some  q*  G  Q. 

Note  that  the  observations  in  (32)  are  defined  to  simulate  measurements  from  a 
single  individual.  The  experimental  data.  [38],  however,  include  measurements 
of  TCE  concentrations  in  several  different  rats.  To  approximate  this  inter- 
individual  variability  within  observations,  we  generate  two  additional  types 
of  simulated  data..  The  first  type  of  data  (Type  I)  represents  the  case  where 
measurements  are  collected  from  several  individuals  over  time,  so  that  each  in¬ 
dividual  is  measured  at  each  time  point.  Note  that  this  type  of  inter-individual 
data  includes  trajectories  over  time  for  each  individual  in  the  group. 

The  second  type  of  data  (Type  11)  represents  the  case  where  measurements 
are  collected  from  multiple  individuals  so  that  each  individual  is  utilized  only 
once.  That  is,  at  each  time  point  there  is  a  separate  group  of  individuals 
that  is  measured.  This  is  the  type  of  data  that  emerges  from  experimental 
measurements  of  tissue  concentrations  when  the  animal  must  be  sacrificed 
and  the  entire  tissue  is  removed  for  analysis.  The  experimental  data  collected 
by  Evans  et  al.  [38]  are  of  this  type. 

For  each  of  these  types  of  data,  we  generate  simulated  observations  by  assum¬ 
ing  that  the  model  parameters  vary  across  the  population  so  that  the  param- 
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eter  vector  q  has  a  probability  distribution  function  P  G  V(Q).  We  assume 
there  are  a  total  of  Ns  individuals  sampled  at  each  time  point  i  =  1 , ,Nt. 
Thus  for  Type  I  data  there  is  a  total  of  Ns  individuals,  while  for  Type  II  data 
there  is  a  total  of  NsNt  individuals. 

For  Type  I  data  we  define  the  observations 


N, 


y  X  r  € 
JVs  3=1 


(33) 


for  i  =  1,  where  q*  is  sampled  from  a  given  P*  G  V(Q)  for  j  = 

1, ...  ,NS.  These  observations  represent  the  expected  value  of  measurements 
taken  over  the  Ns  individuals  in  the  group,  and  are  realizations  of  the  obser¬ 
vations  used  in  (3)  and  described  in  (4). 


The  Type  II  observations  are  given  by 


n. 


(34) 


s  l 


for  i  =  1 , ...  ,Nt,  where  q*j  is  sampled  from  P*  for  j  =  1 , ...  ,NS  and  %  = 


Nt. 


> 


4  -2. 2  Parameter  estimation  results 

We  have  conducted  a  thorough  study  in  [36]  of  the  deterministic  parameter 
estimation  problem  (30)  for  the  adipose  model  parameters 


Q—  [Db,  Dj,  Da,  Hj a,  Hba,  Ubi,  ri,  AB  ■ 


(35) 


We  considered  estimation  problems  for  single  parameters  as  well  as  pairs, 
triples,  and  the  entire  8-vector  q.  The  data  used  in  the  estimation  problems 
included  simulated  data  approximating  a  single  individual  (32)  and  Type  II 
data  (34).  The  objective  function  given  in  (30)  was  minimized  in  Matlab  with 
fminsearch,  which  uses  a  Nelder-Mead  direct  search  algorithm.  A  complete 
description  of  the  estimation  problems  and  results  is  presented  in  [36].  As  the 
probabilistic  estimation  methods  are  the  focus  of  this  paper,  here  we  merely 
summarize  the  results  from  the  deterministic  estimation  method  for  compar¬ 
ative  purposes. 

Using  observations  (32)  simulating  a  single  individual,  we  obtained  optimized 
parameters  that  were  significantly  close  to  the  data-generating  parameters 
q*,  with  more  accurate  estimates  for  pairs  and  triples  of  parameters  than  for 
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the  entire  eight-dimensional  vector  q.  We  studied  the  effect  of  varying  the 
time  points  ti,i  =  1, . . . ,  iV*  at  which  the  observations  are  “collected”  and 
the  number  of  samples  Ns  taken  at  each  time  point,  determining  that  the 
accuracy  of  our  parameter  estimates  increased  as  the  number  of  time  points 
Nt  increased,  as  the  range  of  time  included  in  the  observations  increased,  and 
as  the  number  of  samples  Ns  at  each  time  point  increased. 

When  we  introduced  inter-individual  variability  into  the  observations  by  us¬ 
ing  Type  II  data  (34),  the  deterministic  estimation  method  yielded  parameters 
that  were  less  accurate  as  the  degree  of  population  variability  increased.  These 
results,  which  are  presented  in  full  detail  in  [36] ,  suggested  that  the  determin¬ 
istic  parameter  estimation  technique  may  not  be  the  best  method  to  use  when 
the  observations  are  subject  to  significant  inter- individual  variability. 

To  further  illustrate  this  point,  here  we  present  results  for  the  determinis¬ 
tic  parameter  estimation  problem  with  observations  that  simulate  parameters 
with  a  bimodal  probability  distribution.  In  this  case  we  utilize  the  estimation 
problem  (30)  with  q  =  T>b,  the  dispersion  coefficient.  This  parameter  is  a  mea¬ 
sure  of  the  degree  of  variability  that  occurs  within  an  individual’s  fat  tissue, 
and  it  is  plausible  to  assume  that  T>B  may  be  bimodally  distributed  over  a 
population  with  male  and  female  subpopulations. 

The  observations  used  in  these  estimation  problems  were  generated  with  a 
bimodal  distribution  P*  =  Pf„  composed  of  two  normal  distributions  with 
means  /i\  =  1,  qt2  =  3,  standard  deviations  oq  =  0.1667,  oq  =  0.2  and  mixing 
parameter  0.5  (i.e. ,  equal  weighting  between  the  two  gaussians).  See  Figure  3 
for  a  graph  of  the  probability  density  functions  that  are  combined  to  create  the 
bimodal  density  function.  This  probability  density  is  used  to  generate  obser¬ 
vations  of  Type  I  (33)  and  Type  II  (34);  single  individual  observations  (32)  are 
generated  using  q*  —  1.  We  used  three  different  time  vectors  tt , i  =  1 , ,Nt 
for  our  simulated  observations  to  study  the  effect  of  the  quantity  of  data  on  the 
quality  of  the  estimation  results.  The  vectors  we  used  included  the  following: 


fi  =  [0,  5,  20,  40,  60,  120] 

t2  =  [ 0,  5,  20,  40,  60,  120,  180,  240,  300] 

t3  =  [0,  5,  20,  40,  60,  90,  105,  110,  115,  120,  125,  130, 

135,  150,  180,  210,  240,  270,  275,  285,  290,  295,  300] 

(in  minutes)  with  tf  —  2  hours,  Nt  —  6  for  t3,  tf  =  5  hours,  Nt  —  9  for  t2,  and 
tf  —  5  hours,  Nt  =  23  for  t3.  Simulated  observations  (32)  -  (34)  were  then 
generated  using  these  time  vectors.  The  initial  iterate  go  used  in  the  optimizer 
was  sampled  from  a  uniform  distribution  on  [0.75,1.25]  and  has  a  value  of 
go  =  0.9577. 
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Bimodal  probability  density  function  with  gaussian  components 
2.5 1 - 1 - 1 - 1 - 1 - 1 - 1 - 


Fig.  3.  Normal  probability  density  functions  with  means  Hi  =  1,  fi2  =  3  and  stan¬ 
dard  deviations  <J\  =  0.1667,  <72  =  0.2. 


R 

£2 

£3 

Qopt 

1.0000 

1.0000 

1.0000 

J{qo,z) 

0.0966 

0.5689 

2.2229 

J ( Qopti  z) 

2.5014e— 9 

5.7127e  — 9 

1.4811e  — 8 

Table  1 

Optimized  solutions  qopt  for  the  deterministic  estimation  problem  (30)  with  initial 
iterate  qo  =  0.9577  and  data-generating  parameter  q*  =  1.  The  observations  z  were 
generated  to  simulate  a  single  individual  as  in  (32)  at  the  time  points  t\.  £2  and  1 3. 

For  observations  simulating  a  single  individual,  the  optimized  solutions  were 
equal  to  qopt  =  1.0000  for  each  choice  of  time  vectors  £1,  £2,  £3.  See  Table  1  for 
a  summary  of  these  parameter  estimation  results. 

We  also  carried  out  the  deterministic  estimation  problem  using  Type  I  data  as 
in  (33)  with  the  time  vectors  £1 ,  £2  and  £3,  where  q*  is  sampled  from  the  bimodal 
distribution  function  P*  —  for  j  —  1, . . . ,  Ns.  The  values  of  the  number  Ns 
of  samples  taken  at  each  time  point  we  considered  included  Ns  =  5, 10,  20  and 
50.  Results  are  presented  in  Table  2,  and  plots  of  the  resulting  concentration 
profiles  are  given  in  Figures  4-7.  I11  this  and  all  following  figures,  the  small 
solid  dots  are  the  individual  data  points  and  the  larger  open  circles  are  the 
averaged  values  of  the  data  at  each  point  in  time.  Note  that  the  data  used  in 
the  optimization  problem  are  the  larger  open  circles. 

As  seen  in  the  figures,  as  the  value  of  Ns  increases,  the  predicted  TCE  adipocyte 
concentrations  for  qopt  more  closely  approximate  the  response  from  the  data- 
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u 

t2 

h 

Ns 

5 

10 

20 

50 

5 

10 

20 

50 

5 

10 

20 

50 

Qopt 

2.82 

2.38 

1.79 

1.60 

2.82 

2.32 

1.71 

1.53 

2.82 

2.33 

1.72 

1.54 

J{qo,z) 

34.5 

27.0 

15.1 

10.9 

145 

114 

62.9 

45.5 

593 

464 

257 

186 

J {Qopt-)  z) 

8e  — 6 

4e— 3 

0.02 

0.02 

2e— 4 

0.02 

0.10 

0.12 

8e— 4 

0.08 

0.45 

0.53 

Table  2 

Optimized  solutions  qopt  for  the  deterministic  estimation  problem  (30)  with  initial 
iterate  qo  =  0.9577  and  data-generating  parameter  q*  =  1.  The  observations  2  were 
generated  to  simulate  Type  I  data  as  in  (33)  at  the  time  points  t\.  t?  and  t$  with 
Ns  =  5,10,20  and  50. 

generating  parameters  q*.  For  the  case  Ns  —  5,  note  that  the  individual  data 
points  all  are  clustered  near  the  response  for  the  parameter  value  T>b  =  3, 
so  that  observations  for  this  particular  estimation  problem  do  not  behave 
as  bimodally  distributed  data.  This  is  reflected  in  the  optimized  parameter 
qopt  =  2.82,  which  is  significantly  different  from  the  results  for  other  data 
sets  which  contain  points  from  both  gaussians.  Increasing  the  number  of  time 
points  at  which  the  observations  are  collected  does  not  improve  the  results  for 
this  type  of  data  (e.g.,  see  Figure  7).  These  results  indicate  the  importance  of 
ensuring  a  large  enough  sample  size  for  the  observations. 

Similar  results  for  Type  II  data  are  given  in  Table  3  and  Figures  10  -  13. 
Note  that  as  Ns  increases  and  as  the  number  Nt  of  time  points  increases,  the 
deterministic  method  yields  optimized  parameters  that  generate  reasonable 
approximations  for  the  dynamics  of  TCE  as  seen  with  the  data-generating 
parameters.  It  is  clear  from  Tables  2  and  3,  however,  that  the  optimized  pa¬ 
rameters  qopt  do  not  provide  any  information  about  the  population  variance 
or  the  distribution  of  parameter  values.  Indeed,  most  of  the  values  of  qopt  lie 
in  the  range  from  1.3  to  1.8,  which  is  in  between  the  means  of  the  two  gaus¬ 
sians  that  form  the  bimodal  distribution  P*  =  P^.  Moreover,  as  constants, 
the  optimized  parameters  cannot  suggest  the  underlying  bimodal  nature  of 
the  observations  used  in  the  optimization  problem.  Therefore  it  appears  that 
the  deterministic  approach  may  not  be  most  appropriate  except  in  situations 
when  the  form  of  the  parameter  distribution  function  is  already  known  and  it 
is  necessary  only  to  estimate  the  mean. 


f.3  Probabilistic  parameter  estimation  methods  for  the  TCE  model 


In  this  section  we  apply  the  probability-based  parameter  estimation  methods 
presented  in  Section  2  to  the  TCE  PBPK-hybrid  model.  These  results  are  then 
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h 

t-2 

h 

Ns 

5 

10 

20 

50 

5 

10 

20 

50 

5 

10 

20 

50 

Qopt 

1.34 

1.44 

1.33 

1.28 

1.51 

1.49 

1.52 

1.46 

1.35 

1.42 

1.50 

1.50 

J{qo,z ) 

5.50 

7.26 

5.23 

4.11 

46.5 

42.7 

46.6 

39.7 

128 

148 

174 

174 

J ( Qopt’)  z) 

0.52 

0.12 

0.34 

0.14 

2.49 

0.88 

1.82 

0.98 

21.6 

13.2 

6.08 

3.71 

Table  3 

Optimized  solutions  qopt  for  the  deterministic  estimation  problem  (30)  with  initial 
iterate  qo  =  0.9577  and  data-generating  parameter  q*  =  1.  The  observations  2  were 
generated  to  simulate  Type  II  data  as  in  (34)  at  the  time  points  t\ ,  t,2  and  t%  with 
Ns  =  5,10,20  and  50. 

compared  to  the  results  for  the  deterministic  estimation  problem  as  given  in 
Section  4.2.  Here  we  utilize  the  estimation  problem  with  objective  function 

Nt 

J(P,z)  =  Y,\£[Oy(ti,q)\P\  -  zt? 

i= 1 


and  its  finite-dimensional  approximation 


Nt 

j(pm,z)=y : 


i= 1 


X  Oyfc  qj)pj  -  * 

3=0 


(36) 


as  discussed  in  Section  2,  with  observations  ij,  i  =  1, . . . ,  Nt  corresponding  to 
the  expected  value  (4). 

As  in  Section  4.2  with  the  deterministic  problem,  we  estimate  the  dispersion 
coefficient  T>b  using  observations  from  a  single  individual  and  from  bimodally 
distributed  data  of  Types  I  and  II.  The  admissible  parameter  space  Q  is  defined 
as  the  closed  interval  [0,4],  and  we  use  the  finite  dimensional  approximation 
spaces  Qm  C  Q  for  M  =  32,  64, 128  with 

q.j  =  4 j/M,  j  =  0, . . . ,  M. 


The  constrained  quadratic  programming  problem  (36)  was  solved  in  Matlab 
using  quadprog  to  obtain  the  probabilities  [p0,  ■  ■  ■  ,Pm\  subject  to  the  con¬ 
straints  pj  >  0  and  YJjLo  P]  =  1-  That  is,  we  minimized 

pTAp  +  2bTp  +  c 

subject  to  pj  >  0,  j  —  0, . . . ,  M  and  J2jLoPj  =  1,  where 
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Nt 

Ajk  =  Y  <ij)Oy{U\  Qk),  j,k  =  o, . . . ,  m 

i= 1 

Nt 

bj  =  ~Y  Oy{U ;  qj)zi,  j  —  0, M 

i— 1 
Nt 

c  =  E42 

with  observations  zt  from  (32),  (33)  or  (34),  respectively,  and  time  vectors 
ti,  t-2  or  f3.  As  in  Section  4.2.1,  the  probability  distribution  P*  we  used  to 
generate  the  observations  is  the  bimodal  distribution  Pf„  with  means  /q  —  1 
and  p2  =  3,  standard  deviations  ay  =  0.1667  and  a2  =  0.3333  and  mixing 
parameter  0.5. 


For  the  estimation  problem  with  observations  (32)  simulating  a  single  indi¬ 
vidual,  the  solution  vector  popt  obtained  by  the  optimization  routine  is  given 

by 


p’opt  = 


1  if 

0  otherwise 


Qj  =  1 


for  j  =  0, . . . ,  M  and  M  =  32,  64, 128.  This  solution  corresponds  exactly 
to  the  data-generating  parameter  distribution  with  a  single  mass  at  T>b  —  1, 
suggesting  that  the  probability-based  estimation  method  is  equally  accurate  in 
solving  for  parameters  from  a  single  individual  as  the  deterministic  estimation 
method.  That  is,  if  the  data  correspond  to  a  deterministic  parameter  system, 
using  the  more  general  probability-based  formulation  will  still  provide  correct 
results  by  returning  a  Dirac  measure. 


Results  for  the  probabilistic  method  using  Type  I  data  are  also  plotted  in 
Figures  4-9.  Each  of  these  figures  illustrates  the  case  with  M  =  32;  the 
results  for  M  =  64  and  M  =  128  were  qualitatively  similar  and  are  not 
presented  here.  Simulated  TCE  adipocyte  concentrations  corresponding  to  the 
optimized  parameters  qprob  =  Popt  are  depicted  in  Figures  4  -  7  in  comparison 
with  the  observations  and  the  model  response  corresponding  to  the  optimized 
parameter  qaet  from  the  deterministic  problem  for  the  same  observation  set. 
Graphs  of  the  optimized  probability  distributions  popt  are  given  in  Figures  8 
and  9  for  t\  and  increasing  values  of  Ns. 


Note  in  Figures  4,  5  and  6  that  as  the  sample  size  Ns  increases  from  5  to 
50,  the  predicted  adipocyte  concentrations  for  both  qdet  and  qprob  more  closely 
match  the  response  corresponding  to  the  data-generating  parameter  set.  Sim¬ 
ilarly,  the  optimized  probability  distributions  appear  to  converge  to  the  data- 
generating  distribution  P*  =  Py  as  Ns  increases  (see  Figures  8  and  9).  This 
apparent  convergence  is  in  agreement  with  the  established  theoretical  conver¬ 
gence  of  the  finite-dimensional  parameter  estimation  solutions  summarized  in 
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Section  3  and  detailed  in  [33]. 


As  discussed  in  Section  4.2,  however,  Figure  7  suggests  that  the  accuracy  of 
the  predicted  model  responses  does  not  appear  to  improve  for  this  type  of 
data  as  the  number  Nt  of  time  points  increases  without  a  concurrent  increase 
in  Ns.  In  this  case,  the  additional  time-course  information  contained  in  the 
extra  time  points  does  not  contribute  to  the  richness  of  data  since  the  same 
group  of  individuals  is  sampled  at  every  time  point. 

Unlike  the  optimized  solutions  from  the  deterministic  problem,  the  solutions 
from  the  probabilistic  estimation  method  contain  information  about  the  dis¬ 
tribution  of  the  parameters  across  the  population.  In  addition  to  providing 
an  estimate  of  the  mean  and  variance,  the  probabilistic  method  also  yields  an 
approximation  of  the  probability  distribution  itself.  It  is  important,  however, 
to  ensure  that  a  large  enough  sample  size  is  used  for  the  observations  so  that 
the  population  is  accurately  represented  in  the  data. 

We  present  results  for  the  probabilistic  estimation  problem  with  Type  II  data 
in  Figures  10  -  15.  Predicted  TCE  adipocyte  concentrations  are  given  in  Fig¬ 
ures  10  -  13  for  t\  and  f3  and  various  values  of  Ns ;  probability  distribution 
functions  for  f3  and  varying  Ns  are  presented  in  Figures  14  and  15. 

As  illustrated  in  Figures  10  and  11  for  t\  with  Ns  =  5  and  20  respectively, 
the  optimized  solutions  for  the  probabilistic  method  can  yield  highly  inac¬ 
curate  predictions  of  TCE  adipocyte  concentrations  while  the  deterministic 
method  produces  an  accurate  match  to  the  observations.  Our  studies  with 
the  quadratic  programming  problem  have  indicated  the  existence  of  multiple 
solutions  that  satisfy  the  constraints,  with  some  of  the  solutions  yielding  inac¬ 
curate  predictions  of  adipocyte  concentrations  beyond  the  time  period  covered 
in  the  data. 

This  difficulty  does  not  arise  for  observations  that  utilize  the  larger  time  vector 
f3  (see  Figures  12  and  13),  where  the  predicted  responses  for  both  methods 
are  reasonably  close  to  the  observations  for  all  four  values  of  Ns.  As  seen  in 
Figures  14  and  15,  however,  the  optimized  probability  distributions  Popt  do 
not  appear  to  converge  with  increasing  values  of  Ns  as  clearly  as  in  the  case 
with  Type  I  data. 

It  is  important  to  note  that  the  quality  of  the  Type  II  data  is  significantly  dif¬ 
ferent  than  the  quality  of  Type  I  data  since  there  is  no  time-course  information 
from  any  individual  in  the  Type  II  data.  Since  each  individual  is  measured 
only  once  for  the  Type  II  data,  the  data  points  do  not  contain  any  trajectories 
in  time  from  an  individual.  In  contrast,  every  individual  is  measured  at  each 
time  point  for  the  Type  I  data,  incorporating  important  time-course  dynamics 
into  the  observations. 
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As  demonstrated  in  our  results,  the  deterministic  method  can  reasonably  pre¬ 
dict  the  expected  value  of  the  model  response  over  the  population  for  Type 
I  and  Type  II  data,  but  can  provide  no  information  about  the  population 
distribution  of  the  parameters.  The  probabilistic  method,  however,  can  suc¬ 
cessfully  yield  the  expected  value,  variance  and  overall  shape  of  the  population 
distribution  using  Type  I  data.  The  probabilistic  method  has  some  difficulty 
capturing  this  probability  distribution  with  Type  II  data  since  valuable  time- 
course  information  is  not  contained  in  these  data.  Therefore  it  appears  that 
Type  II  data  is,  not  surprisingly,  less  desirable  for  situations  when  the  prob¬ 
ability  distribution  of  parameters  must  be  estimated  with  no  prior  knowledge 
of  the  shape  of  the  distribution. 


5  Concluding  remarks 


In  this  paper  we  have  presented  probability-based  parameter  estimation  meth¬ 
ods  which  incorporate  and  account  for  uncertainty  and  population  variability 
that  arise  in  biological  models.  We  outlined  known  theoretical  results  that 
address  issues  of  well-posedness  for  these  estimation  problems,  and  we  applied 
them  to  a  toxicokinetic  model  for  trichloroethylene  using  simulated  obser¬ 
vations  that  exhibited  differing  types  of  variability.  These  results  were  com¬ 
pared  to  results  obtained  with  a  traditional  deterministic  parameter  estimation 
method. 

As  one  might  expect,  our  results  indicate  that  the  performance  of  the  deter¬ 
ministic  and  probabilistic  estimation  methods  depends  greatly  on  both  the 
quantity  and  quality  of  the  data  used  in  the  estimation  process.  For  data 
that  represent  a  single  individual,  each  method  produced  parameter  estimates 
that  accurately  matched  the  data-generating  parameter  q*.  Not  surprisingly, 
in  this  case  it  appears  that  the  deterministic  method  is  sufficient  for  estimating 
parameters  using  data  from  a  single  individual. 

When  any  degree  of  variability  is  introduced  into  the  experimental  data,  how¬ 
ever,  the  deterministic  estimation  method  is  unable  to  capture  any  of  this 
variability  in  its  parameter  estimates.  This  method  can  in  some  cases  rea¬ 
sonably  predict  the  expected  value  of  the  parameter  over  the  population,  but 
can  provide  no  information  about  the  variance  or  the  shape  of  the  probability 
distribution. 

For  observations  of  the  expected  value  that  contain  important  time-course 
trajectories  (as  in  the  case  of  our  Type  I  data),  the  probabilistic  method  is 
able  to  successfully  predict  both  the  expected  value  and  the  overall  probability 
distribution  of  the  parameter.  This  is  accomplished  by  solving  a  standard 
quadratic  programming  problem  with  no  prior  knowledge  of  the  population 
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distribution.  As  seen  in  our  results,  it  is  important  to  use  a  sufficiently  large 
sample  size  so  that  the  observations  contain  an  adequate  representation  of  the 
population. 


The  results  for  the  probability-based  estimation  problem  are  less  conclusive 
for  Type  If  data,  which  are  significantly  different  from  the  Type  1  data  in 
that  they  contain  only  one  observation  from  each  individual.  This  results  in 
data  that  include  no  individual  time-course  trajectories,  which  contributes  to 
difficulties  in  the  optimization  process.  There  appear  to  be  many  solutions  to 
the  constrained  quadratic  programming  problem  for  this  type  of  data,  and 
some  of  these  solutions  produce  inaccurate  predictions  of  TCE  adipocyte  con¬ 
centrations  for  time  periods  not  included  in  the  data  themselves.  Increasing 
the  number  of  time  points  Nt  at  which  the  observations  are  collected  seems  to 
improve  the  results,  but  the  convergence  of  the  probability  distributions  is  less 
clear  than  for  the  case  with  Type  I  data.  Our  results  therefore  suggest  that 
it  is  less  than  desirable  to  utilize  Type  II  data  for  estimating  parameters  that 
have  an  unknown  probability  distribution.  If  at  all  possible,  it  is  best  to  col¬ 
lect  experimental  data  that  contain  multiple  measurements  from  individuals 
over  time,  as  this  time-course  information  adds  significant  richness  to  the  ex¬ 
perimental  data  and  significantly  enhances  our  ability  to  estimate  underlying 
inter-individual  variability. 


A  key  feature  of  the  probability-based  methods  presented  here  is  their  ability 
to  estimate  population  distributions  for  parameters  without  the  use  of  priors. 
In  situations  where  there  is  little  information  about  the  shape  of  the  probabil¬ 
ity  distribution(s),  this  can  be  a  clear  advantage  over  methods  (e.g.,  Markov 
Chain-Monte  Carlo)  that  require  the  specification  of  priors.  For  example,  a  pa¬ 
rameter  with  a  bimodal  parameter  distribution  may  be  mistakenly  estimated 
as  a  unimodal  distribution  if  too  much  weight  is  placed  on  the  assumption  of  a 
unimodal  prior.  The  probability-based  estimation  methods  avoid  this  problem, 
and  the  utilization  of  Dirac  delta  measures  in  the  proper  metric  space  guaran¬ 
tees  the  theoretical  convergence  of  the  resulting  estimates  of  the  probability 
distributions. 


Current  and  future  efforts  related  to  this  work  include  a  study  and  comparison 
of  other  probability-based  parameter  estimation  techniques.  In  particular,  we 
plan  to  implement  Markov  Chain-Monte  Carlo  methods  applied  to  a  PBPK 
model  and  compare  the  results  we  obtain  with  the  probability-based  estima¬ 
tions  outlined  in  this  paper. 
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Bimodal,  qM  =  32,  tf  =  2,  Nt  =  6,  Ns  =  5  (Type  I) 


Fig.  4.  Simulated  observations  and  predicted  TCE  adipocyte  concentrations  us¬ 
ing  optimized  parameters  q^et  from  the  deterministic  estimation  problem  (30)  and 
Qprob  =  Popt  from  the  probabilistic  problem  (36)  with  M  =  32.  Observations  used 
in  the  estimation  problems  were  Type  I  data  with  t\  (tf  =  2  hours,  Nt  =  6)  and 
with  Ns  =  5.  In  this  and  all  similar  figures  that  follow,  the  solid  black  dots  are 
individual  observations  and  the  black  open  circles,  which  are  used  in  the  estimation 
problem,  are  the  averaged  values  of  the  individual  observations  at  each  time  point. 
The  solid  line  is  the  model  response  corresponding  to  the  data-generating  bimodal 
distribution  P*  =  using  (4)  with  P* .  The  dashed  and  dotted  lines  are  the  model 
responses  corresponding  to  qdet  and  qprob ,  respectively. 
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Bimodal,  qM  =  32,  tf  =  2,  Nt  =  6,  Ns  =  1 0  (Type  I) 


Fig.  5.  Simulated  observations  and  predicted  TCE  adipocyte  concentrations  using 
optimized  parameters  q^t  from  the  deterministic  estimation  problem  (30)  and  qprob 
from  the  probabilistic  problem  (36)  with  M  =  32.  Observations  used  in  the  estima¬ 
tion  problems  were  Type  I  data  with  t\  (tf  =  2  hours,  N  =  6)  and  with  Ns  =  10. 
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Bimodal,  qM  =  32,  tf  =  2,  Nt  =  6,  Ns  =  20  (Type  I) 


Fig.  6.  Simulated  observations  and  predicted  TCE  adipocyte  concentrations  using 
optimized  parameters  qdet  from  the  deterministic  estimation  problem  (30)  and  qprob 
from  the  probabilistic  problem  (36)  with  M  =  32.  Observations  used  in  the  estima¬ 
tion  problems  were  Type  I  data  with  t\  (tf  =  2  hours,  Nt  =  6)  and  with  Ns  =  20 
(top  figure)  and  Ns  =  50  (bottom  figure). 
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Bimodal,  qM  =  32,  tf  =  5,  Nt  =  9,  Ns  =  1 0  (Type  I) 


Fig.  7.  Simulated  observations  and  predicted  TCE  adipocyte  concentrations  using 
optimized  parameters  qdet  from  the  deterministic  estimation  problem  (30)  and  qprob 
from  the  probabilistic  problem  (36)  with  M  =  32.  Observations  used  in  the  esti¬ 
mation  problems  were  Type  I  data  with  Ns  =  10  and  with  t,2  (top  figure)  and  t% 
(bottom  figure). 
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Bimodal,  qM  =  32,  tf  =  2,  Nt  =  6,  Ns  =  5  (Type  1) 
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Bimodal,  qM  =  32,  tf  =  2,  Nt  =  6,  Ns  =  10  (Type  I) 
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Fig.  8.  Probability  distribution  functions  P  as  a  function  of  the  parameter  q  =  T>b 
for  the  data-generating  distribution  P*  =  P^  and  the  optimized  distribution  Pprob 
from  the  probabilistic  estimation  problem  (36)  with  Type  I  data,  M  =  32,  ti  and 
with  Ns  =  5  (top  figure)  and  Ns  =  10  (bottom  figure). 
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Bimodal,  qM  =  32,  tf  =  2,  Nt  =  6,  Ns  =  20  (Type  I) 
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Bimodal,  qM  =  32,  tf  =  2,  Nt  =  6,  Ns  =  50  (Type  I) 
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Fig.  9.  Probability  distribution  functions  P  as  a  function  of  the  parameter  q  =  T>b 
for  the  data-generating  distribution  P*  =  P^  and  the  optimized  distribution  Pprob 
from  the  probabilistic  estimation  problem  (36)  with  Type  I  data,  M  =  32,  ti  and 
with  Ns  =  20  (top  figure)  and  Ns  =  50  (bottom  figure). 


35 


Bimodal,  qM  =  32,  tf  =  2,  Nt  =  6,  Ns  =  5  (Type  II) 


Fig.  10.  Simulated  observations  and  predicted  TCE  adipocyte  concentrations  using 
optimized  parameters  qdet  from  the  deterministic  estimation  problem  (30)  and  qprob 
from  the  probabilistic  problem  (36)  with  M  =  32.  Observations  used  in  the  estima¬ 
tion  problems  were  Type  II  data  with  t\  (tf  =  2  hours,  Nt  =  6)  and  with  Ns  =  5 
(top  figure)  and  Ns  =  10  (bottom  figure). 
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Bimodal,  qM  =  32,  tf  =  2,  Nt  =  6,  Ns  =  20  (Type  II) 


Fig.  11.  Simulated  observations  and  predicted  TCE  adipocyte  concentrations  using 
optimized  parameters  qdet  from  the  deterministic  estimation  problem  (30)  and  qprob 
from  the  probabilistic  problem  (36)  with  M  =  32.  Observations  used  in  the  estima¬ 
tion  problems  were  Type  II  data  with  t\  {tf  =  2  hours,  Nt  =  6)  and  with  Ns  =  20 
(top  figure)  and  Ns  =  50  (bottom  figure). 
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Bimodal,  qM  =  32,  tf  =  5,  Nt  =  23,  Ns  =  5  (Type  II) 


Fig.  12.  Simulated  observations  and  predicted  TCE  adipocyte  concentrations  using 
optimized  parameters  q^et  from  the  deterministic  estimation  problem  (30)  and  qprob 
from  the  probabilistic  problem  (36)  with  M  =  32.  Observations  used  in  the  estima¬ 
tion  problems  were  Type  II  data  with  t3  [tf  =  5  hours,  Nt  =  23)  and  with  Ns  =  5 
(top  figure)  and  Ns  =  10  (bottom  figure). 
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Bimodal,  qM  =  32,  tf  =  5,  Nt  =  23,  Ns  =  20  (Type  II) 


Fig.  13.  Simulated  observations  and  predicted  TCE  adipocyte  concentrations  using 
optimized  parameters  qdet  from  the  deterministic  estimation  problem  (30)  and  qprob 
from  the  probabilistic  problem  (36)  with  M  =  32.  Observations  used  in  the  estima¬ 
tion  problems  were  Type  II  data  with  1 3  {tf  =  5  hours,  Nt  =  23)  and  with  Ns  =  20 
(top  figure)  and  Ns  =  50  (bottom  figure). 
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Bimodal,  qM  =  32,  tf  =  5,  Nt  =  23,  Ns  =  5  (Type  II) 
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Fig.  14.  Probability  distribution  functions  P  as  a  function  of  the  parameter  q  =  T>b 
for  the  data-generating  distribution  P*  =  P^  and  the  optimized  distribution  Pprob 
from  the  probabilistic  estimation  problem  (36)  with  Type  II  data,  M  =  32,  ts  and 
with  Ns  =  5  (top  figure)  and  Ns  =  10  (bottom  figure). 
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Bimodal,  qM  =  32,  tf  =  5,  Nt  =  23,  Ns  =  20  (Type  II) 
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Fig.  15.  Probability  distribution  functions  P  as  a  function  of  the  parameter  q  =  T>b 
for  the  data-generating  distribution  P*  =  P^  and  the  optimized  distribution  Pprob 
from  the  probabilistic  estimation  problem  (36)  with  Type  II  data,  M  =  32,  ts  and 
with  N 5  =  20  (top  figure)  and  Ns  =  50  (bottom  figure). 
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